rtsVis Overview
Johannes_Mast
12 11 2020
rtsVis
rtsVis is a package with a single purpose: The creation of animated visualizations of time satellite imagery time series.
In this overview, we show how to create such visualizations on the example of the MODIS MCD43A4 product of the Western Cape, South Africa.
Getting Started
To get started, we load a number of packages. To begin with, we need rtsVis itself. We use raster to load and plot rasters and sf for vector objects.
rtsVis also goes hand in hand with moveVis from which we use some core functions. We also use magrittr for the convenient pipes.
library(rtsVis)
library(raster)
library(sf)
library(moveVis)
library(magrittr) Perhaps most importantly, we need the data itself: A time series of raster images.
Preparing the Time Series
rtsVis works on lists of objects - lists of raster objects or lists of frames (plots). Most functions of the package take a list as input, modify it, and return the modified list.
The first step in a rtsVis pipeline is to get our time series into this format.
We can load the images using the raster::stack function.
modis_filenames <- list.files("Data/Modis_images/",pattern = ".tif$",full.names = T)
r_list <- lapply(modis_filenames,stack)Or maybe we have a time series prepared in some other format. Sometimes, you will already have done some preprocessing, such as resampling or masking. What matters is that the images are matching in projection and extent. We have such a time series already stored as .rds file.
r_list <- readRDS("Beispieldaten/MODIS/WesternCape/MODIS_WesternCape.rds")First, we take a look at one of our images.
plotRGB(r_list[[20]],1,4,3,stretch="lin")
We can see that it looks quite good. But there are very small patches of what appear to be missing values.
Filling NAs
To fix them, we use the first rtsVis function: ts_fill_na(). This uses the raster::approxNA() function to fill NA values from the closest temporal neighbor in the time series - a crude approach, but reasonably quick.
r_list_filled <- rtsVis::ts_fill_na(r_list)
plotRGB(r_list_filled[[20]],1,4,3,stretch="lin")
Much better! In the southwest, out on the Atlantic Ocean, we have an area which cannot be filled because we have no values for it in any of our MODIS images.
Interpolating the time series
Our MODIS Time Series is already at a very good temporal resolution of one day. Great for analysis - but for visualizations, we do not want to be restricted by the temporal resolution of our source.
We will use the rtsVis::ts_raster function to linearly interpolate our time series to a higher temporal resolution. For this, we need:
- List of input images
r_list - acquisition times
r_timesfor each input image - times
out_timesfor each desired output image
Our input images we have already. So all we need are the times.
Input times
How we get the times depends on the source of our images. Here, we stored our acquisition dates in the names of our list object. In other cases, we might read them out of the metadata or the file names. However we get there, what matters is that we get a vector of POSIX objects of the same length as our list of images, containing the input time for each image.
head(names(r_list_filled))## [1] "2019-07-01" "2019-07-02" "2019-07-03" "2019-07-04" "2019-07-05"
## [6] "2019-07-06"in_dates <- as.POSIXct(names(r_list_filled))
head(in_dates)## [1] "2019-07-01 CEST" "2019-07-02 CEST" "2019-07-03 CEST" "2019-07-04 CEST"
## [5] "2019-07-05 CEST" "2019-07-06 CEST"Output times
Our output times are all the dates for which we will receive an interpolated output image. The output times, therefore, should also be a vector of POSIX (although they need not match the input images). How we get them is, again, up to us. If we only want to increase the temporal resolution, we can simply create a sequence of POSIX, like this:
out_dates <-seq.POSIXt(from = in_dates[1],
to = in_dates[length(in_dates)],
length.out = length(in_dates)*2 )
head(out_dates)## [1] "2019-07-01 00:00:00 CEST" "2019-07-01 11:57:51 CEST"
## [3] "2019-07-01 23:55:42 CEST" "2019-07-02 11:53:33 CEST"
## [5] "2019-07-02 23:51:24 CEST" "2019-07-03 11:49:15 CEST"Interpolating
We pass our images and times to ts_raster. The option fade_raster activates linear interpolation instead of nearest-neighbor interpolation.
r_list_filled_interpolated <- rtsVis::ts_raster(r_list = r_list_filled,
r_times = in_dates,
out_times = out_dates,
fade_raster = T)length(r_list_filled)## [1] 154length(r_list_filled_interpolated)## [1] 308We have doubled the length of our time series. We should bear in mind that the interpolated values no longer correspond to actual measurements and should not be used for analysis. Also, any further processing will take longer. So we should use this power responsibly. For visualizations, however, this is a great boon, as we have doubled the frequency of our images.
Also, the resulting list has the time written in its attributes. This is crucial for the application of some rtsVis functions.
head(rtsVis:::.ts_get_frametimes(r_list_filled_interpolated))## [1] "2019-07-01 00:00:00 CEST" "2019-07-01 11:57:51 CEST"
## [3] "2019-07-01 23:55:42 CEST" "2019-07-02 11:53:33 CEST"
## [5] "2019-07-02 23:51:24 CEST" "2019-07-03 11:49:15 CEST"Creating the visualisations
Now we have a time series in our required format and at our desired resolution. We can go on to create some visualization.
All the visualizations we create will be ggplots. From our list of raster objects will come a list of ggplot frames, one for each input raster object.
- Flow frames, visualizing the distribution of raster values at certain positions such as points or polygons in the form of graphs. Different types of graphs are available:
- violin plots
- line plots
- custom plots
- Spatial frames, created by
ts_makeframes(). These can be further enhanced by adding- Geographical style objects such as north arrows, using
moveVis::add_northarrow() - Spatial positions such as points or polygons using
rtsVis::ts_add_positions_to_frames()
- Geographical style objects such as north arrows, using
The returned list of frames can be animated using moveVis::animate_frames().
Flow Frames
Flow frames visualize amounts, distributions or proportions of raster values within an area. They can, in principle, take many forms such as lines, histograms, boxplots, scatterplots, pie charts or barcharts. rtsVis currently implements violin plots and line plots. However, custom plot types can be used as well.
By default, these plots will be calculated for the entire area covered by the raster. However, if we want to distinguish between different areas of interest within our study area. In this case, we can provide the spatial objects and the plots will use the spatial objects as a grouping variable.
Violin Plots
Here, we use four administrative areas which lie within out study area.
polygons <- st_read("Beispieldaten/Ancillary/WesternCape/Western_Cape_polygons.shp") %>% st_transform(crs = st_crs(r_list_filled_interpolated[[1]]))
#Plot the polygons
plotRGB(r_list_filled_interpolated[[50]],1,4,3,stretch="lin")
plot(st_geometry(polygons),add=T,col="red",alpha=0.5)
We can see that they are of different sizes and shapes and cover different land uses. This may be important when evaluating our plots later. But first, we create them. Only the first three arguments are truly necessary, however, we make use of the option to customize the plots.
violin_frames <- rtsVis::ts_flow_frames(r_list = r_list_filled_interpolated,
positions = polygons,
plot_function = "violin",
# Optional Arguments
position_names = polygons$NAME_3,
band_names = c("620 - 670","841 - 876","459 - 479","545 - 565"),
band_colors = c("firebrick3","darkorchid3","dodgerblue3","olivedrab3"),
band_legend_title = "Wavelength [nm]",
position_legend = F,
legend_position = "left",
band_legend=F,
aes_by_pos = F)The result is a list of the same length as our list of input rasters.
length(violin_frames)## [1] 308And each object in our new list is a ggplot. We can take a look at an individual violin_frame.
violin_frames[[50]]## Warning: Removed 68 rows containing non-finite values (stat_ydensity).
Line Plots
In the same way we used ts_flow_frames() to create a violin plot of polygon data, we can use it to create a line plot of data under a set of points.
points <- st_read("Beispieldaten/Ancillary/WesternCape/Western_Cape_pts.shp") %>% st_transform(crs = st_crs(r_list_filled_interpolated[[1]]))
line_frames <- rtsVis::ts_flow_frames(r_list = r_list_filled_interpolated,
positions = points[c(1,3,4),],
plot_function = "line",
# Optional Arguments
position_names = points$pointname[c(1,3,4)],
band_names = c("620 - 670","841 - 876","459 - 479","545 - 565"),
band_colors = c("firebrick3","darkorchid3","dodgerblue3","olivedrab3"),
band_legend_title = "Wavelength [nm]",
position_legend = T,
legend_position = "left",
band_legend=T,
aes_by_pos = T)
line_frames[[length(line_frames)]]## Reading layer `Western_Cape_pts' from data source `D:\DEV\TSVis\Beispieldaten\Ancillary\WesternCape\Western_Cape_pts.shp' using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 2 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 1682063 ymin: -3736239 xmax: 1797136 ymax: -3663320
## projected CRS: Sinusoidal
We could, in principle, also create a violin plot of the data extracted under points. Then we would probably want to use a buffer to our point data, like so:
points <- st_read("Beispieldaten/Ancillary/WesternCape/Western_Cape_pts.shp") %>% st_transform(crs = st_crs(r_list_filled_interpolated[[1]]))
line_frames <- rtsVis::ts_flow_frames(r_list = r_list_filled_interpolated,
positions = points[c(1,3,4),],
plot_function = "violin",
pbuffer = 1337, ##Use values in a circular buffer of 1337
##map units around the points
#
position_names = points$pointname[c(1,3,4)],
band_names = c("620 - 670","841 - 876","459 - 479","545 - 565"),
band_colors = c("firebrick3","darkorchid3","dodgerblue3","olivedrab3"),
band_legend_title = "Wavelength [nm]",
position_legend = T,
legend_position = "left",
band_legend=F,
aes_by_pos = T)
line_frames[[length(line_frames)]]Custom Plot Types
Violin and line charts are nice for some applications but often other plot types may be more suitable. Besides the options “line” and “violin”, we can also provide ts_flow_frames() with a custom plotting function.
There are a few requirements for such a function. It should:
- usually subset the extracted tibble dataframe to be in line with the temporal scope
- return a ggplot object
- accept (and optionally, but ideally also implement) the following arguments, which are prepared by ts_flow_frames:
- i: The index of the current frame. This is used to grab the data which is relevant for the current frame from edf
- edf: The extracted tibble containing values for plotting. The format of this tibble will be discussed in detail in the following section.
- pl: Position legend. A logical value that indicates whether a legend for the positions should be plotted. Only useful if abp that is, if any aesthetic is mapped to position.
- bl: Band legend. A logical value that indicates whether a legend for the bands should be plotted.
- plt: Position legend title.
- blt: Band legend title.
- lp: Legend position.
- ps: Plot size. This should be passed to a ggplot geom. As it is treated differently depending on the type of geom, it should sometimes be adjusted by a factor.
- vs: Value sequence. A sequence of values like (0,1250,2500,3750,5000). This is calculated across all extracted values and allows the user to pass breaks in a value scale.
- abp: Aesthetics by position? A logical value that indicates whether some aesthetic should be mapped to the position. For example, every point could be represented by a different linetype in a line plot.
Below, we create a plotting function that creates a barchart.
library(ggplot2)
barplot <- function(i, edf , pl,lp, bl, blt,plt, ps, vs,abp){
#The data up to the current frame (this will be plotted)
x = edf[edf$frame == i,]
#All data (this sets the frame)
y=edf
## generate base plot
p <- ggplot(x, aes(x = band,y=value,fill=band))
## style it
p <- p +
geom_bar(stat = "identity")+
theme_bw() +
theme(aspect.ratio = 1) +
scale_y_continuous( breaks = vs,limits = c(min(vs),max(vs)))+
theme(legend.position = lp)+
facet_grid(~position_name)+
scale_x_discrete(guide = guide_axis(n.dodge = 2))
#add the colors
p <- p +
scale_fill_manual(values = x$band_colors,breaks = x$band, name=blt)
return(p)
}We can then pass this function to ts_flow_frames().
bar_frames <- rtsVis::ts_flow_frames(r_list = r_list_filled_interpolated,
positions = points[c(1,3,4),],
plot_function = barplot,
position_names = points$pointname[c(1,3,4)],
band_names = c("620 - 670","841 - 876","459 - 479","545 - 565"),
band_colors = c("firebrick3","darkorchid3","dodgerblue3","olivedrab3"),
band_legend_title = "Wavelength [nm]",
position_legend = T,
legend_position = "left",
band_legend=T,
aes_by_pos = T)
bar_frames[[50]]
Many other plot types are possible. For inspiration, visit the r-graph-gallery!
Also visit Claus Wilke’s excellent book which deals with the most common pitfalls of data visualization.
Spatial Frames
Spatial frames are, quite simply, images that can be discrete, gradient or RGB, depending of the type of input raster.
Creating the Spatial Frames
The frames can be created using the ts_makeframes function.
Two things to note here:
- Since we are working with MODIS images, we specify the bands we need to get a true color composite. If we were jusing a different data source, or wanted a false color composite, we could adjust our band selection accordingly.
- We apply a a quantile-based stretch. The quantiles will be derived from values across the entire time series and therefore the stretch will be consistent.
r_frames <- rtsVis::ts_makeframes(x_list = r_list_filled_interpolated,
l_indices = c(1,4,3),
minq = 0.2,
maxq = 0.999)## [1] "Guessing raster type."
## [1] "Detected 3+ layers, choosing raster type 'RGB'"r_frames[[50]] # Examine a single frame
We can see that our frames are plotted as ggplots with latitude and longitude as the axes.
Adding additional Elements to the Frames
We can use functionality from the moveVis package to enhance our images with additional elements.
r_frames_styled <- r_frames %>%
moveVis::add_labels(x = "Longitude", y = "Latitude")%>%
moveVis::add_northarrow(colour = "white", position = "bottomright") %>%
moveVis::add_timestamps(type = "label") %>%
moveVis::add_progress()
r_frames_styled[[50]] # Examine a single frame
Now the spatial and temporal context of the map is greatly increased. In some cases, we want to go even further.
Adding Positions to the Frames
We may want to visualize positions, certain points or areas of significance on the map. This can be useful to visualize the areas of interest we explore in the chats. In our case here, we use the same positions we previously used to create the flow frames.
In our first case, we add the municipalities - which we used to create the violin plots - and add them as polygons to our frames.
r_frames_with_polygons <-
rtsVis::ts_add_positions_to_frames(r_frame_list = r_frames_styled,
positions = polygons,
aes_by_pos = T,
legend_position = "right",
position_names = polygons$NAME_3,
position_legend_title = "Municipality",
psize = 1,
add_text = T,
tsize = 3,
t_hjust = 5000,
ttype = "label")
r_frames_with_polygons[[50]] # Examine a single frame
In our second case, we do it with our previously used points.
r_frames_with_points <-
rtsVis::ts_add_positions_to_frames(r_frame_list = r_frames_styled,
position_names = points$pointname,
positions = points,
add_text = T,
tsize = 4,
ttype = "label",
t_vjust = 8500)
r_frames_with_points[[50]] # Examine a single frame
Joining Frames
we can use moveVis::join_frames to join frames. This way we can combine our charts which show the distribution of values within our areas of interest with a spatial element that presents their spatial context. See the join_frames() documentation for more details.
point_joined<- moveVis::join_frames(frames_lists = list(r_frames_with_points,line_frames))
point_joined[[10]]
Animating Frames
Lists of frames, as returned by ts_makeframes(), ts_flow_frames(), or moveVis::join_frames(), can be animated in a number of different formats ùsingmoveVis function, moveVis::animate_frames(). See the documentation of the function for more details. The resulting file is written to disk.
moveVis::animate_frames(point_joined,
out_file = "WesternCape_MODIS_point.gif",
width=1600)